import datashader as ds
import datashader.transfer_functions as tf
import datashader.glyphs
from datashader import reductions
from datashader.core import bypixel
from datashader.utils import lnglat_to_meters as webm, export_image
from datashader.colors import colormap_select, Greys9, viridis, inferno
import copy
from pyproj import Proj, transform
import numpy as np
import pandas as pd
import urllib
import json
import datetime
import colorlover as cl
import plotly.offline as py
import plotly.graph_objs as go
from plotly import tools
from shapely.geometry import Point, Polygon, shape
# In order to get shapley, you'll need to run [pip install shapely.geometry] from your terminal
from functools import partial
from IPython.display import GeoJSON
py.init_notebook_mode()
For module 2 we'll be looking at techniques for dealing with big data. In particular binning strategies and the datashader library (which possibly proves we'll never need to bin large data for visualization ever again.)
To demonstrate these concepts we'll be looking at the PLUTO dataset put out by New York City's department of city planning. PLUTO contains data about every tax lot in New York City.
PLUTO data can be downloaded from here. Unzip them to the same directory as this notebook, and you should be able to read them in using this (or very similar) code. Also take note of the data dictionary, it'll come in handy for this assignment.
# Code to read in v17, column names have been updated (without upper case letters) for v18
# bk = pd.read_csv('PLUTO17v1.1/BK2017V11.csv')
# bx = pd.read_csv('PLUTO17v1.1/BX2017V11.csv')
# mn = pd.read_csv('PLUTO17v1.1/MN2017V11.csv')
# qn = pd.read_csv('PLUTO17v1.1/QN2017V11.csv')
# si = pd.read_csv('PLUTO17v1.1/SI2017V11.csv')
# ny = pd.concat([bk, bx, mn, qn, si], ignore_index=True)
ny = pd.read_csv('nyc_pluto_18v2_csv/pluto_18v2.csv', low_memory=False) # Warning said to use 'low_memory=False' but I really have no idea why...
# Getting rid of some outliers
ny = ny[(ny['yearbuilt'] > 1850) & (ny['yearbuilt'] < 2020) & (ny['numfloors'] != 0)]
I'll also do some prep for the geographic component of this data, which we'll be relying on for datashader.
You're not required to know how I'm retrieving the lattitude and longitude here, but for those interested: this dataset uses a flat x-y projection (assuming for a small enough area that the world is flat for easier calculations), and this needs to be projected back to traditional lattitude and longitude.
wgs84 = Proj("+proj=longlat +ellps=GRS80 +datum=NAD83 +no_defs")
nyli = Proj("+proj=lcc +lat_1=40.66666666666666 +lat_2=41.03333333333333 +lat_0=40.16666666666666 +lon_0=-74 +x_0=300000 +y_0=0 +ellps=GRS80 +datum=NAD83 +to_meter=0.3048006096012192 +no_defs")
ny['xcoord'] = 0.3048*ny['xcoord']
ny['ycoord'] = 0.3048*ny['ycoord']
ny['lon'], ny['lat'] = transform(nyli, wgs84, ny['xcoord'].values, ny['ycoord'].values)
ny = ny[(ny['lon'] < -60) & (ny['lon'] > -100) & (ny['lat'] < 60) & (ny['lat'] > 20)]
#Defining some helper functions for DataShader
background = "black"
export = partial(export_image, background = background, export_path="export")
cm = partial(colormap_select, reverse=(background!="black"))
Part 1: Binning and Aggregation¶
Binning is a common strategy for visualizing large datasets. Binning is inherent to a few types of visualizations, such as histograms and 2D histograms (also check out their close relatives: 2D density plots and the more general form: heatmaps.
While these visualization types explicitly include binning, any type of visualization used with aggregated data can be looked at in the same way. For example, lets say we wanted to look at building construction over time. This would be best viewed as a line graph, but we can still think of our results as being binned by year:
trace = go.Scatter(
# I'm choosing BBL here because I know it's a unique key.
x = ny.groupby('yearbuilt').count()['bbl'].index,
y = ny.groupby('yearbuilt').count()['bbl']
)
layout = go.Layout(
xaxis = dict(title = 'Year Built'),
yaxis = dict(title = 'Number of Lots Built')
)
fig = go.Figure(data = [trace], layout = layout)
py.iplot(fig)
Something looks off... You're going to have to deal with this imperfect data to answer this first question.
But first: some notes on pandas. Pandas dataframes are a different beast than R dataframes, here are some tips to help you get up to speed:
Hello all, here are some pandas tips to help you guys through this homework:
Indexing and Selecting: .loc and .iloc are the analogs for base R subsetting, or filter() in dplyr
Group By: This is the pandas analog to group_by() and the appended function the analog to summarize(). Try out a few examples of this, and display the results in Jupyter. Take note of what's happening to the indexes, you'll notice that they'll become hierarchical. I personally find this more of a burden than a help, and this sort of hierarchical indexing leads to a fundamentally different experience compared to R dataframes. Once you perform an aggregation, try running the resulting hierarchical datafrome through a reset_index().
Reset_index: I personally find the hierarchical indexes more of a burden than a help, and this sort of hierarchical indexing leads to a fundamentally different experience compared to R dataframes. reset_index() is a way of restoring a dataframe to a flatter index style. Grouping is where you'll notice it the most, but it's also useful when you filter data, and in a few other split-apply-combine workflows. With pandas indexes are more meaningful, so use this if you start getting unexpected results.
Indexes are more important in Pandas than in R. If you delve deeper into the using python for data science, you'll begin to see the benefits in many places (despite the personal gripes I highlighted above.) One place these indexes come in handy is with time series data. The pandas docs have a huge section on datetime indexing. In particular, check out resample, which provides time series specific aggregation.
Merging, joining, and concatenation: There's some overlap between these different types of merges, so use this as your guide. Concat is a single function that replaces cbind and rbind in R, and the results are driven by the indexes. Read through these examples to get a feel on how these are performed, but you will have to manage your indexes when you're using these functions. Merges are fairly similar to merges in R, similarly mapping to SQL joins.
Apply: This is explained in the "group by" section linked above. These are your analogs to the plyr library in R. Take note of the lambda syntax used here, these are anonymous functions in python. Rather than predefining a custom function, you can just define it inline using lambda.
Browse through the other sections for some other specifics, in particular reshaping and categorical data (pandas' answer to factors.) Pandas can take a while to get used to, but it is a pretty strong framework that makes more advanced functions easier once you get used to it. Rolling functions for example follow logically from the apply workflow (and led to the best google results ever when I first tried to find this out and googled "pandas rolling")
Google Wes Mckinney's book "Python for Data Analysis," which is a cookbook style intro to pandas. It's an O'Reilly book that should be pretty available out there.
Question¶
After a few building collapses, the City of New York is going to begin investigating older buildings for safety. The city is particularly worried about buildings that were unusually tall when they were built, since best-practices for safety hadn’t yet been determined. Create a graph that shows how many buildings of a certain number of floors were built in each year (note: you may want to use a log scale for the number of buildings). Find a strategy to bin buildings (It should be clear 20-29-story buildings, 30-39-story buildings, and 40-49-story buildings were first built in large numbers, but does it make sense to continue in this way as you get taller?)
Binning Data¶
I spent an insane amount of time trying different methods of binning and grouping data and had a lot of trouble getting the data into the format needed in order to be able to plot in plotly. I left a lot of the code below even though some of it wasn't really needed in order to show what I tried and keep as notes for future reference.
I binned data for the number of floors and for the decade in which a building was constructed and added those as new columns to the original dataframe as well as creating grouped aggragated dataframes with the same bins.
I discoverd after MANY failures that in order to plot the data in plotly any categorical variables MUST have labels. Unlabeled columns just would not plot even if coverted to a list or series. Indexes also had to be converted to columns with labels if they were categorical.
Number of Floors¶
# Start your answer here, inserting more cells as you go along
#ny[['numfloors','yearbuilt']][1:5]
#ny.groupby('yearbuilt').count()['bbl'].index
bins1 = [0, 1, 2, 3, 4, 5, 10, 15, 20, 25, 30, 35, 40, 45, 60, 75, 150, 300]
labels1 = ['1', '2', '3', '4', '5', '6-10', '11-15', '16-20', '21-25', '26-30', '31-35', '36-40', '41-45', '46-60', '61-75', '76-150', '151-300']
floors = pd.cut(ny['numfloors'], bins1, labels=labels1)
ny['numfloors_bins'] = floors
ny['numfloors_bins'].head()
Changed the code above and then updated the chart below to test different bin sizes and find appropriate cutoffs¶
Floors = pd.DataFrame(ny['numfloors'].groupby(floors).count())
Floors.columns = ['count']
Floors.reset_index(inplace=True)
Floors
Test plot to make sure this works...¶
I also wanted to get an idea of the distribution of different size buildings. The plot makes it clear that it will probably be necessary to log transform the data to see the distributions later since the bin for 2 story buildings is more than twice the size of any of the others and it can't be cut any smaller.
trace = go.Bar(
# I'm choosing BBL here because I know it's a unique key.
x = ny.groupby('numfloors_bins').count()['bbl'].index,
y = ny.groupby('numfloors_bins').count()['bbl']
)
layout = go.Layout(
xaxis = dict(title = 'numfloors_bins'),
yaxis = dict(title = 'Number of Lots Built')
)
fig = go.Figure(data = [trace], layout = layout)
py.iplot(fig)
Year Built¶
bins2 = range(1850, 2021, 10)
labels2 = []
for year in range(1850, 2020, 10):
labels2.append(str(year) + "'s")
dates = pd.cut(ny['yearbuilt'], bins2, labels=labels2)
ny['yearbuilt_bins'] = pd.cut(ny['yearbuilt'], bins2, labels=labels2)
ny['yearbuilt_bins'].head()
Changed the code above and then updated the chart below to test different bin sizes and find appropriate cutoffs¶
Year_Built = pd.DataFrame(ny['yearbuilt'].groupby(dates).count())
Year_Built.columns = ['count']
Year_Built.reset_index(inplace=True)
Year_Built
Another test plot to make sure this works...¶
You can see a huge building boom in the 1920's!
trace = go.Bar(
# I'm choosing BBL here because I know it's a unique key.
x = ny.groupby('yearbuilt_bins').count()['bbl'].index,
y = ny.groupby('yearbuilt_bins').count()['bbl']
)
layout = go.Layout(
xaxis = dict(title = 'Year Built'),
yaxis = dict(title = 'Number of Lots Built')
)
fig = go.Figure(data = [trace], layout = layout)
py.iplot(fig)
Check the new columns in our dataframe...¶
ny.iloc[:, -2:].head(10)
Create new aggragated dataframe¶
After trying unsuccessfully many times to use the new categorical columns that I added to the original ny dataframe, I gave up and created a new dataframe with just the aggragated data for plotting... Even this gave me problems until I finally figured out how to reset the index here: https://github.com/pandas-dev/pandas/issues/19136
data = pd.DataFrame(ny['bbl'].groupby([floors, dates]).count().unstack())
#data.reset_index(inplace=True) Throws error!!!!
# Got the function below here: https://github.com/pandas-dev/pandas/issues/19136
def reset_index(df):
'''Returns DataFrame with index as columns'''
index_df = df.index.to_frame(index=False)
df = df.reset_index(drop=True)
# In merge is important the order in which you pass the dataframes
# if the index contains a Categorical.
# pd.merge(df, index_df, left_index=True, right_index=True) does not work
return pd.merge(index_df, df, left_index=True, right_index=True)
data = reset_index(data)
data
Plot data for each year¶
After an insane amopunt of time trying all different kinds of plots and trying to tweak the labels etc... This is the best I could come up with so far. I still can't figure out how to fix a bunch of things that are bugging me:
- Should have axis labels on the bottom and down the let side of the entire grid, but when IO tried to add them they wound up just on the first plot. (I removed them since they were not working right)
- Can't figure out how to leave the last spot in the grid blank. It automatically repeated the first trace there even though I didn't assign anything to that spot.
- I'm not loving the log transfromaion, but it is harder to see bars for the taller buildings with the non-log transformed data
- Can't figure out how to get ALL the plots to have the same y-axis, not just the same on each row... As a result it's difficult to see the building boom in the 20's. All the plots look like they have similar maximum numbers if you aren't paying attention to the axes.
- I also would have liked to have the plots go three across, but when I had 3 columns every other label on the x-axis was missing because I guess the labels would have run into each other that way and I couldn't figure out how to change the size of the text for the labels.
fig = tools.make_subplots(cols = 2, rows = 9, print_grid = False, shared_yaxes=True,
subplot_titles=(data.columns[1:]))
traces = []
for y in labels2:
traces.append({'type' : 'bar',
'name' : y,
'x' : data['numfloors'].tolist(),
'y' : np.log10(data[y]).tolist()
})
traces.append(traces[0])
t=0
while t < 16:
for r in range(1,10):
for c in range(1,3):
fig.append_trace(trace = traces[t], row = r, col = c)
t +=1
fig['layout'].update(height=1500, width=800, title='Number of Floors - Log10 Transformed Building Counts',
#xaxis = dict(title = 'Number of Floors'), yaxis = dict(title = 'Number of Lots Built'),
showlegend = False)
#fig = go.Figure(data = [traces[2]], layout = layout)
py.iplot(fig)
#help(tools.make_subplots)
fig
Part 2: Datashader¶
Datashader is a library from Anaconda that does away with the need for binning data. It takes in all of your datapoints, and based on the canvas and range returns a pixel-by-pixel calculations to come up with the best representation of the data. In short, this completely eliminates the need for binning your data.
As an example, lets continue with our question above and look at a 2D histogram of YearBuilt vs NumFloors:
yearbins = 200
floorbins = 200
yearBuiltCut = pd.cut(ny['yearbuilt'], np.linspace(ny['yearbuilt'].min(), ny['yearbuilt'].max(), yearbins))
numFloorsCut = pd.cut(ny['numfloors'], np.logspace(1, np.log(ny['numfloors'].max()), floorbins))
xlabels = np.floor(np.linspace(ny['yearbuilt'].min(), ny['yearbuilt'].max(), yearbins))
ylabels = np.floor(np.logspace(1, np.log(ny['numfloors'].max()), floorbins))
data = [
go.Heatmap(z = ny.groupby([numFloorsCut, yearBuiltCut])['bbl'].count().unstack().fillna(0).values,
colorscale = 'Greens', x = xlabels, y = ylabels)
]
py.iplot(data)
This shows us the distribution, but it's subject to some biases discussed in the Anaconda notebook Plotting Perils.
Here is what the same plot would look like in datashader:
cvs = ds.Canvas(800, 500, x_range = (ny['yearbuilt'].min(), ny['yearbuilt'].max()),
y_range = (ny['numfloors'].min(), ny['numfloors'].max()))
agg = cvs.points(ny, 'yearbuilt', 'numfloors')
view = tf.shade(agg, cmap = cm(Greys9), how='log')
export(tf.spread(view, px=2), 'yearvsnumfloors')
That's technically just a scatterplot, but the points are smartly placed and colored to mimic what one gets in a heatmap. Based on the pixel size, it will either display individual points, or will color the points of denser regions.
Datashader really shines when looking at geographic information. Here are the latitudes and longitudes of our dataset plotted out, giving us a map of the city colored by density of structures:
NewYorkCity = (( -74.29, -73.69), (40.49, 40.92))
cvs = ds.Canvas(700, 700, *NewYorkCity)
agg = cvs.points(ny, 'lon', 'lat')
view = tf.shade(agg, cmap = cm(inferno), how='log')
export(tf.spread(view, px=2), 'firery')
Interestingly, since we're looking at structures, the large buildings of Manhattan show up as less dense on the map. The densest areas measured by number of lots would be single or multi family townhomes.
Unfortunately, Datashader doesn't have the best documentation. Browse through the examples from their github repo. I would focus on the visualization pipeline and the US Census Example for the question below. Feel free to use my samples as templates as well when you work on this problem.
Question¶
You work for a real estate developer and are researching underbuilt areas of the city. After looking in the Pluto data dictionary, you've discovered that all tax assessments consist of two parts: The assessment of the land and assessment of the structure. You reason that there should be a correlation between these two values: more valuable land will have more valuable structures on them (more valuable in this case refers not just to a mansion vs a bungalow, but an apartment tower vs a single family home). Deviations from the norm could represent underbuilt or overbuilt areas of the city. You also recently read a really cool blog post about bivariate choropleth maps, and think the technique could be used for this problem.
Datashader is really cool, but it's not that great at labeling your visualization. Don't worry about providing a legend, but provide a quick explanation as to which areas of the city are overbuilt, which areas are underbuilt, and which areas are built in a way that's properly correlated with their land value.
# assuming that this is a valid calculation?
# Makes sense that the total value minus the land value should equal the structure value, right?
ny['assessstructure'] = ny['assesstot'] - ny['assessland']
ny[['assessland', 'assessstructure', 'assesstot']].head(12)
# Bin the data in three equal sized categories for both land value and structure value
val_labels = ['low', 'med', 'high']
ny['land_group'] = pd.qcut(ny['assessland'], 3, labels=val_labels)
ny['structure_group'] = pd.qcut(ny['assessstructure'], 3, labels=val_labels)
ny[['land_group', 'structure_group']].head(12)
ny['land_group'].value_counts()
ny['structure_group'].value_counts()
# Combine the two binned categories into one 9 category variable
ny['combined_value'] = ny['land_group'].str.cat(ny['structure_group'].values.astype(str), sep='-')
ny['combined_value'] = pd.Categorical(ny['combined_value'])
ny['combined_value'].head(12)
# assign colors to the 9 categories
# land values range from gray to bright pink and structure values range from gray to bright blue
colors = {'low-low': '#e8e8e8',
'low-med': '#ace4e4',
'low-high': '#5ac8c8',
'med-low': '#dfb0d6',
'med-med': '#a5add3',
'med-high': '#5698b9',
'high-low': '#be64ac',
'high-med': '#8c62aa',
'high-high': '#3b4994'
}
# Plot the data
NewYorkCity = (( -74.29, -73.69), (40.49, 40.92))
cvs = ds.Canvas(900, 900, *NewYorkCity)
agg = cvs.points(ny, 'lon', 'lat', ds.count_cat('combined_value'))
view = tf.shade(agg, color_key = colors)
export(tf.spread(view, px=2), 'pretty_colors')
Conclusion¶
Almost everything in Manhattan seems to fall within the high-high (dark-blue) category, while the outer boroughs show more variance. There are quite a few areas in the low-med, med-high, and low-high (teal-blue) categories indicating some overbuilt areas in Brookyln and Queens, with far less of them scattered through the Bronx and Staten Island. The North Shore of Staten Island has the largest area in the brightest pink low-high category indicating underbuilt land, but there are also very large areas in the other pink/purple (med-low and high-med) colors that would indicate less extremely underbuilt areas in parts of Brooklyn and Queens that are furthest from Manhattan.
Trying a different approach...¶
Not sure about my method for cutting the data above. If there are extreme outlier values that would actually get hidden by this approach, because I cut the data into basically 3 equal parts for low, med, and high. So anything that is far above or far below normal would just get labeled as high or low. So I thought maybe I should try cutting the data into thirdds based on the max and min values instead...
# Get max values
land_max = ny['assessland'].max()
struc_max = ny['assessstructure'].max()
# Create bins based on max values
land_bins = [0, land_max/3, (land_max/3)*2, land_max]
land_bins
struc_bins = [0, struc_max/3, (struc_max/3)*2, struc_max]
struc_bins
I also checked the minimum and there were a significant number of zeros as well as very low values, so I thought to leave zero as a minimum. The provblem is that the maximum value is so far above the rest of the values that by binning this way almost everything falls in the low category for both land assessment and structure value. You can see the few outliers in the tables below.
ny['land_group2'] = pd.cut(ny['assessland'], bins=land_bins, labels=val_labels)
ny['structure_group2'] = pd.cut(ny['assessstructure'], bins=struc_bins, labels=val_labels)
ny[['assessland', 'assessstructure', 'structure_group2']][ny['land_group2']=='med']
ny[['assessland', 'assessstructure', 'structure_group2']][ny['land_group2']=='high']
ny[['assessland', 'assessstructure', 'structure_group2']][ny['structure_group2']=='med']
ny[['assessland', 'assessstructure', 'structure_group2']][ny['structure_group2']=='high']
Trying a log transformation...¶
So next I thought maybe a log transformation on the assessment varibles first might solve the problem...
# Get log transfromed data, fill NA's and neg infinity with zeros
log_land = np.log(ny['assessland']).fillna(0)
log_structure = np.log(ny['assessstructure']).fillna(0)
from numpy import isneginf
log_land[isneginf(log_land)]=0
log_structure[isneginf(log_structure)]=0
Visualization of the log transformed land values¶
Even after log transforming the data, there is still a strong skew in the land values.
trace = go.Histogram(
x = log_land,
nbinsx = 100
)
layout = go.Layout(
xaxis = dict(title = 'log assessland')
)
fig = go.Figure(data = [trace], layout = layout)
py.iplot(fig)